library("dplyr")
library("ggplot2")
library('gridExtra')
Sys.setlocale("LC_TIME")
[1] "en_US.UTF-8"
df <- read.csv(file = './data/datos_ruido_madrid.csv', header = FALSE, sep = ';', strip.white = TRUE)
dfs <- read.csv(file = './data/EstacionesMedidaControlAcustico.csv', header = FALSE, sep = ';', strip.white = TRUE, fileEncoding = "ISO-8859-1")
colnames(df) <- c("Station","Year","Month","Day","Period","LAeq","LAS01","LAS10","LAS50","LAS90","LAS99")
df$Date <- as.Date(with(df, paste(Year, Month, Day,sep="-")), "%Y-%m-%d")
df$MonthName <- months(df$Date)
df$DayMonthName <- with(df, paste(weekdays(df$Date), MonthName,sep="-"))
df$LAS01 <- NULL
df$LAS10 <- NULL
df$LAS50 <- NULL
df$LAS90 <- NULL
df$LAS99 <- NULL
head(df)
dfs <- read.csv(file = './data/EstacionesMedidaControlAcustico.csv', header = TRUE, sep = ',', strip.white = TRUE, fileEncoding = "ISO-8859-1")
dfs[ ,c('COD_VIA', 'VIA_CLASE','VIA_PAR','VIA_NOMBRE','Dirección','Longitud_gms','Latitud_gms','LATITUD_ED50','LONGITUD_ED50','Alt..m.','Fecha.alta','Coordenada_X_ETRS89','Coordenada_Y_ETRS89','LONGITUD_WGS84','LATITUD_WGS84','X.2','X.3')] <- list(NULL)
dfs <- dfs %>% 
  rename(
    Station=Nº,
    StationName = Nombre,
    Long = X, 
    Lat = X.1
    )
head(dfs)
df <- merge(df, dfs, by = "Station", all.x = TRUE)
head(df)

Calculamos la media de los índices por año, mes, estación y periodo, y añadimos las variables al dataframe.

meanLAeqByMonthStationAndPeriod <- df %>% 
    group_by(Year,MonthName, StationName, Period) %>% 
    summarise(meanLAeqByMonth = mean(LAeq))
df <- merge(df, meanLAeqByMonthStationAndPeriod, by = c("Year","MonthName", "StationName", "Period"), all.x = TRUE)
colnames(df)
 [1] "Year"            "MonthName"       "StationName"     "Period"          "Station"        
 [6] "Month"           "Day"             "LAeq"            "Date"            "DayMonthName"   
[11] "Long"            "Lat"             "meanLAeqByMonth"
head(df)

Nos centramos únicamente en el índice LAeq durante el periodo ‘T’. Nos interesa conocer la media mensual por año de ese índice para cada una de las estaciones de Madrid.

 dfPeriodT <- df %>% filter(Period == "T")
dfMeanLAeqByMonth <- distinct(dfPeriodT, Year, MonthName, StationName, meanLAeqByMonth)
# Year as factor in order to facilitate the choice of colors in the plot
dfMeanLAeqByMonth$Year <- factor(dfMeanLAeqByMonth$Year)
ggplot(dfMeanLAeqByMonth, aes(x = ordered(MonthName,levels=month.name), y = meanLAeqByMonth, group=Year)) +
  geom_point(aes(color=Year), size=2)  +
  geom_line(aes(color=Year))  +
  facet_wrap(~ StationName, ncol = 3, scales="free") +
  labs(title = "Comparación del índice LAeq medido en las distintas estaciones a lo largo de los años",
       subtitle = "Se ha realizado la media mensual de los datos recopilados en cada una de las estaciones",
       y = "Media mensual LAeq") +
 theme(legend.position="top")

Una vez que hemos visto cómo el efecto de la cuarentena debido al coronavirus se refleja en los niveles de ruido, nos centramos únicamente en el periodo de tiempo durante el cual fueron aplicadas las medidas de confinamiento e hibernación de la economía. Para ello estudiaremos los niveles de ruido con un nivel de descripción semanal.

Para ello crearemos obtendremos los promedios diarios de los índices de los años 2015-2019 teniendo en cuenta las siguientes premisas:

  1. Los promedios se tienen que realizar no por fecha sino por días (i.e. se promedian los días de la semana de un mismo mes). Así pues, tomaremos, por ejemplo, todos los lunes de un mismo mes de los distintos años para realizar el promedio y obtener un lunes “típico” del mes en cuestión.

  2. Construiremos un mes “típico” para compararlo con el mes correspondiente de 2020.

  3. El efecto de Madrid central observado en la gráfica anterior no va a ser tenido en cuenta ya que vamos a promediar sobre todos los años previos al de la crisis del coronavirus.

head(df)

Calculamos la media por día del mes para todos los años comprendidos entre 2015 y 2019 (ambos incluidos).

 df_2015_2019<- df %>% filter(Year!=2020)
 meanLAeqByDayMonthStationAndPeriod_2015_2019 <- df_2015_2019 %>%
  group_by(DayMonthName, StationName, Period) %>% 
  summarise(meanLAeqByDayMonth_2015_2019 = mean(LAeq))
head(meanLAeqByDayMonthStationAndPeriod_2015_2019)

Creamos el dataframe para el año 2020 y añadimos el resultado de la media del índice LAeq promediada por tipo de día del mes. Cortamos el dataframe para quedarnos únicamente con los días comprendidos entre los meses marzo-mayo, ambos incluidos. Además, solo nos interesa el Periodo ‘T’.

df2020 <- df %>% filter(Year==2020) %>% filter(Period=='T')
dfCovidEffect <- merge(df2020, meanLAeqByDayMonthStationAndPeriod_2015_2019, by= c("DayMonthName", "StationName", "Period"), x.ALL = True)
dfCovidEffect <- dfCovidEffect %>% filter(Month>2, Month<6) 
dfCovidEffect <- dfCovidEffect[, c('Date', 'StationName','LAeq', 'meanLAeqByDayMonth_2015_2019')]
head(dfCovidEffect)
dfCovidEffectMelted <- reshape2::melt(dfCovidEffect, id.var=c('Date', 'StationName'))
ggplot(dfCovidEffectMelted, aes(x = Date, y = value, col=variable)) +
  geom_point(color='black',size=1)  +
  geom_line()  +
  facet_wrap(~ StationName, ncol = 3, scales="free") +
  labs(title = "Comparación del índice de ruido LAeq diario para 2020 y la media de los años 2015-2019",
       subtitle = "Subtítulo", y = "LAeq") +
 theme(legend.position="top")

LS0tCnRpdGxlOiAiQW7DoWxpc2lzIGRlbCBydWlkbyBlbiBNYWRyaWQiCm91dHB1dDogCiAgICBodG1sX25vdGVib29rOiBkZWZhdWx0IAogICAgaHRtbF9kb2N1bWVudDogZGVmYXVsdAoKLS0tCmBgYHtyfQpsaWJyYXJ5KCJkcGx5ciIpCmxpYnJhcnkoImdncGxvdDIiKQpsaWJyYXJ5KCdncmlkRXh0cmEnKQpgYGAKCmBgYHtyfQpTeXMuc2V0bG9jYWxlKCJMQ19USU1FIikKCmBgYAoKCmBgYHtyfQpkZiA8LSByZWFkLmNzdihmaWxlID0gJy4vZGF0YS9kYXRvc19ydWlkb19tYWRyaWQuY3N2JywgaGVhZGVyID0gRkFMU0UsIHNlcCA9ICc7Jywgc3RyaXAud2hpdGUgPSBUUlVFKQpkZnMgPC0gcmVhZC5jc3YoZmlsZSA9ICcuL2RhdGEvRXN0YWNpb25lc01lZGlkYUNvbnRyb2xBY3VzdGljby5jc3YnLCBoZWFkZXIgPSBGQUxTRSwgc2VwID0gJzsnLCBzdHJpcC53aGl0ZSA9IFRSVUUsIGZpbGVFbmNvZGluZyA9ICJJU08tODg1OS0xIikKCmNvbG5hbWVzKGRmKSA8LSBjKCJTdGF0aW9uIiwiWWVhciIsIk1vbnRoIiwiRGF5IiwiUGVyaW9kIiwiTEFlcSIsIkxBUzAxIiwiTEFTMTAiLCJMQVM1MCIsIkxBUzkwIiwiTEFTOTkiKQoKZGYkRGF0ZSA8LSBhcy5EYXRlKHdpdGgoZGYsIHBhc3RlKFllYXIsIE1vbnRoLCBEYXksc2VwPSItIikpLCAiJVktJW0tJWQiKQoKZGYkTW9udGhOYW1lIDwtIG1vbnRocyhkZiREYXRlKQoKZGYkRGF5TW9udGhOYW1lIDwtIHdpdGgoZGYsIHBhc3RlKHdlZWtkYXlzKGRmJERhdGUpLCBNb250aE5hbWUsc2VwPSItIikpCgpkZiRMQVMwMSA8LSBOVUxMCmRmJExBUzEwIDwtIE5VTEwKZGYkTEFTNTAgPC0gTlVMTApkZiRMQVM5MCA8LSBOVUxMCmRmJExBUzk5IDwtIE5VTEwKCgpoZWFkKGRmKQpgYGAKCmBgYHtyfQoKZGZzIDwtIHJlYWQuY3N2KGZpbGUgPSAnLi9kYXRhL0VzdGFjaW9uZXNNZWRpZGFDb250cm9sQWN1c3RpY28uY3N2JywgaGVhZGVyID0gVFJVRSwgc2VwID0gJywnLCBzdHJpcC53aGl0ZSA9IFRSVUUsIGZpbGVFbmNvZGluZyA9ICJJU08tODg1OS0xIikKCmRmc1sgLGMoJ0NPRF9WSUEnLCAnVklBX0NMQVNFJywnVklBX1BBUicsJ1ZJQV9OT01CUkUnLCdEaXJlY2Npw7NuJywnTG9uZ2l0dWRfZ21zJywnTGF0aXR1ZF9nbXMnLCdMQVRJVFVEX0VENTAnLCdMT05HSVRVRF9FRDUwJywnQWx0Li5tLicsJ0ZlY2hhLmFsdGEnLCdDb29yZGVuYWRhX1hfRVRSUzg5JywnQ29vcmRlbmFkYV9ZX0VUUlM4OScsJ0xPTkdJVFVEX1dHUzg0JywnTEFUSVRVRF9XR1M4NCcsJ1guMicsJ1guMycpXSA8LSBsaXN0KE5VTEwpCgpkZnMgPC0gZGZzICU+JSAKICByZW5hbWUoCiAgICBTdGF0aW9uPU7CuiwKICAgIFN0YXRpb25OYW1lID0gTm9tYnJlLAogICAgTG9uZyA9IFgsIAogICAgTGF0ID0gWC4xCiAgICApCgpoZWFkKGRmcykKCmBgYAoKCmBgYHtyfQpkZiA8LSBtZXJnZShkZiwgZGZzLCBieSA9ICJTdGF0aW9uIiwgYWxsLnggPSBUUlVFKQpoZWFkKGRmKQoKYGBgCgoKQ2FsY3VsYW1vcyBsYSBtZWRpYSBkZSBsb3Mgw61uZGljZXMgcG9yIGHDsW8sIG1lcywgZXN0YWNpw7NuIHkgcGVyaW9kbywgeSBhw7FhZGltb3MgbGFzIHZhcmlhYmxlcyBhbCBkYXRhZnJhbWUuCgpgYGB7cn0KCm1lYW5MQWVxQnlNb250aFN0YXRpb25BbmRQZXJpb2QgPC0gZGYgJT4lIAogICAgZ3JvdXBfYnkoWWVhcixNb250aE5hbWUsIFN0YXRpb25OYW1lLCBQZXJpb2QpICU+JSAKICAgIHN1bW1hcmlzZShtZWFuTEFlcUJ5TW9udGggPSBtZWFuKExBZXEpKQoKZGYgPC0gbWVyZ2UoZGYsIG1lYW5MQWVxQnlNb250aFN0YXRpb25BbmRQZXJpb2QsIGJ5ID0gYygiWWVhciIsIk1vbnRoTmFtZSIsICJTdGF0aW9uTmFtZSIsICJQZXJpb2QiKSwgYWxsLnggPSBUUlVFKQoKY29sbmFtZXMoZGYpCgpoZWFkKGRmKQpgYGAKCk5vcyBjZW50cmFtb3Mgw7puaWNhbWVudGUgZW4gZWwgw61uZGljZSBMQWVxIGR1cmFudGUgZWwgcGVyaW9kbyAnVCcuIE5vcyBpbnRlcmVzYSBjb25vY2VyIGxhIG1lZGlhIG1lbnN1YWwgcG9yIGHDsW8gZGUgZXNlIMOtbmRpY2UgcGFyYSBjYWRhIHVuYSBkZSBsYXMgZXN0YWNpb25lcyBkZSBNYWRyaWQuIAoKYGBge3J9CiBkZlBlcmlvZFQgPC0gZGYgJT4lIGZpbHRlcihQZXJpb2QgPT0gIlQiKQoKZGZNZWFuTEFlcUJ5TW9udGggPC0gZGlzdGluY3QoZGZQZXJpb2RULCBZZWFyLCBNb250aE5hbWUsIFN0YXRpb25OYW1lLCBtZWFuTEFlcUJ5TW9udGgpCgojIFllYXIgYXMgZmFjdG9yIGluIG9yZGVyIHRvIGZhY2lsaXRhdGUgdGhlIGNob2ljZSBvZiBjb2xvcnMgaW4gdGhlIHBsb3QKZGZNZWFuTEFlcUJ5TW9udGgkWWVhciA8LSBmYWN0b3IoZGZNZWFuTEFlcUJ5TW9udGgkWWVhcikKCmBgYAoKCgoKYGBge3IsIGZpZy53aWR0aD0xMiwgZmlnLmhlaWdodD04fQoKZ2dwbG90KGRmTWVhbkxBZXFCeU1vbnRoLCBhZXMoeCA9IG9yZGVyZWQoTW9udGhOYW1lLGxldmVscz1tb250aC5uYW1lKSwgeSA9IG1lYW5MQWVxQnlNb250aCwgZ3JvdXA9WWVhcikpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvcj1ZZWFyKSwgc2l6ZT0yKSAgKwogIGdlb21fbGluZShhZXMoY29sb3I9WWVhcikpICArCiAgZmFjZXRfd3JhcCh+IFN0YXRpb25OYW1lLCBuY29sID0gMywgc2NhbGVzPSJmcmVlIikgKwogIGxhYnModGl0bGUgPSAiQ29tcGFyYWNpw7NuIGRlbCDDrW5kaWNlIExBZXEgbWVkaWRvIGVuIGxhcyBkaXN0aW50YXMgZXN0YWNpb25lcyBhIGxvIGxhcmdvIGRlIGxvcyBhw7FvcyIsCiAgICAgICBzdWJ0aXRsZSA9ICJTZSBoYSByZWFsaXphZG8gbGEgbWVkaWEgbWVuc3VhbCBkZSBsb3MgZGF0b3MgcmVjb3BpbGFkb3MgZW4gY2FkYSB1bmEgZGUgbGFzIGVzdGFjaW9uZXMiLAogICAgICAgeSA9ICJNZWRpYSBtZW5zdWFsIExBZXEiKSArCiB0aGVtZShsZWdlbmQucG9zaXRpb249InRvcCIpCmBgYAoKVW5hIHZleiBxdWUgaGVtb3MgdmlzdG8gY8OzbW8gZWwgZWZlY3RvIGRlIGxhIGN1YXJlbnRlbmEgZGViaWRvIGFsIGNvcm9uYXZpcnVzIHNlIHJlZmxlamEgZW4gbG9zIG5pdmVsZXMgZGUgcnVpZG8sIG5vcyBjZW50cmFtb3Mgw7puaWNhbWVudGUgZW4gZWwgcGVyaW9kbyBkZSB0aWVtcG8gZHVyYW50ZSBlbCBjdWFsIGZ1ZXJvbiBhcGxpY2FkYXMgbGFzIG1lZGlkYXMgZGUgY29uZmluYW1pZW50byBlIGhpYmVybmFjacOzbiBkZSBsYSBlY29ub23DrWEuIFBhcmEgZWxsbyBlc3R1ZGlhcmVtb3MgbG9zIG5pdmVsZXMgZGUgcnVpZG8gY29uIHVuIG5pdmVsIGRlIGRlc2NyaXBjacOzbiBzZW1hbmFsLiAKClBhcmEgZWxsbyBjcmVhcmVtb3Mgb2J0ZW5kcmVtb3MgbG9zIHByb21lZGlvcyBkaWFyaW9zIGRlIGxvcyDDrW5kaWNlcyBkZSBsb3MgYcOxb3MgMjAxNS0yMDE5IHRlbmllbmRvIGVuIGN1ZW50YSBsYXMgc2lndWllbnRlcyBwcmVtaXNhczogCgoxLiBMb3MgcHJvbWVkaW9zIHNlIHRpZW5lbiBxdWUgcmVhbGl6YXIgbm8gcG9yIGZlY2hhIHNpbm8gcG9yIGTDrWFzIChpLmUuIHNlIHByb21lZGlhbiBsb3MgZMOtYXMgZGUgbGEgc2VtYW5hIGRlIHVuIG1pc21vIG1lcykuIEFzw60gcHVlcywgdG9tYXJlbW9zLCBwb3IgZWplbXBsbywgdG9kb3MgbG9zIGx1bmVzIGRlIHVuIG1pc21vIG1lcyBkZSBsb3MgZGlzdGludG9zIGHDsW9zIHBhcmEgcmVhbGl6YXIgZWwgcHJvbWVkaW8geSBvYnRlbmVyIHVuIGx1bmVzICJ0w61waWNvIiBkZWwgbWVzIGVuIGN1ZXN0acOzbi4gCgoyLiBDb25zdHJ1aXJlbW9zIHVuIG1lcyAidMOtcGljbyIgcGFyYSBjb21wYXJhcmxvIGNvbiBlbCBtZXMgY29ycmVzcG9uZGllbnRlIGRlIDIwMjAuIAoKMy4gRWwgZWZlY3RvIGRlIE1hZHJpZCBjZW50cmFsIG9ic2VydmFkbyBlbiBsYSBncsOhZmljYSBhbnRlcmlvciBubyB2YSBhIHNlciB0ZW5pZG8gZW4gY3VlbnRhIHlhIHF1ZSB2YW1vcyBhIHByb21lZGlhciBzb2JyZSB0b2RvcyBsb3MgYcOxb3MgcHJldmlvcyBhbCBkZSBsYSBjcmlzaXMgZGVsIGNvcm9uYXZpcnVzLiAKYGBge3J9CmhlYWQoZGYpCmBgYAoKQ2FsY3VsYW1vcyBsYSBtZWRpYSBwb3IgZMOtYSBkZWwgbWVzIHBhcmEgdG9kb3MgbG9zIGHDsW9zIGNvbXByZW5kaWRvcyBlbnRyZSAyMDE1IHkgMjAxOSAoYW1ib3MgaW5jbHVpZG9zKS4KCmBgYHtyfQoKIGRmXzIwMTVfMjAxOTwtIGRmICU+JSBmaWx0ZXIoWWVhciE9MjAyMCkKCiBtZWFuTEFlcUJ5RGF5TW9udGhTdGF0aW9uQW5kUGVyaW9kXzIwMTVfMjAxOSA8LSBkZl8yMDE1XzIwMTkgJT4lCiAgZ3JvdXBfYnkoRGF5TW9udGhOYW1lLCBTdGF0aW9uTmFtZSwgUGVyaW9kKSAlPiUgCiAgc3VtbWFyaXNlKG1lYW5MQWVxQnlEYXlNb250aF8yMDE1XzIwMTkgPSBtZWFuKExBZXEpKQoKaGVhZChtZWFuTEFlcUJ5RGF5TW9udGhTdGF0aW9uQW5kUGVyaW9kXzIwMTVfMjAxOSkKYGBgCgpDcmVhbW9zIGVsIGRhdGFmcmFtZSBwYXJhIGVsIGHDsW8gMjAyMCB5IGHDsWFkaW1vcyBlbCByZXN1bHRhZG8gZGUgbGEgbWVkaWEgZGVsIMOtbmRpY2UgTEFlcSBwcm9tZWRpYWRhIHBvciB0aXBvIGRlIGTDrWEgZGVsIG1lcy4gQ29ydGFtb3MgZWwgZGF0YWZyYW1lIHBhcmEgcXVlZGFybm9zIMO6bmljYW1lbnRlIGNvbiBsb3MgZMOtYXMgY29tcHJlbmRpZG9zIGVudHJlIGxvcyBtZXNlcyBtYXJ6by1tYXlvLCBhbWJvcyBpbmNsdWlkb3MuIEFkZW3DoXMsIHNvbG8gbm9zIGludGVyZXNhIGVsIFBlcmlvZG8gJ1QnLiAKYGBge3J9CgpkZjIwMjAgPC0gZGYgJT4lIGZpbHRlcihZZWFyPT0yMDIwKSAlPiUgZmlsdGVyKFBlcmlvZD09J1QnKQoKZGZDb3ZpZEVmZmVjdCA8LSBtZXJnZShkZjIwMjAsIG1lYW5MQWVxQnlEYXlNb250aFN0YXRpb25BbmRQZXJpb2RfMjAxNV8yMDE5LCBieT0gYygiRGF5TW9udGhOYW1lIiwgIlN0YXRpb25OYW1lIiwgIlBlcmlvZCIpLCB4LkFMTCA9IFRydWUpCgpkZkNvdmlkRWZmZWN0IDwtIGRmQ292aWRFZmZlY3QgJT4lIGZpbHRlcihNb250aD4yLCBNb250aDw2KSAKCmRmQ292aWRFZmZlY3QgPC0gZGZDb3ZpZEVmZmVjdFssIGMoJ0RhdGUnLCAnU3RhdGlvbk5hbWUnLCdMQWVxJywgJ21lYW5MQWVxQnlEYXlNb250aF8yMDE1XzIwMTknKV0KCmhlYWQoZGZDb3ZpZEVmZmVjdCkKCmBgYAoKYGBge3J9CmRmQ292aWRFZmZlY3RNZWx0ZWQgPC0gcmVzaGFwZTI6Om1lbHQoZGZDb3ZpZEVmZmVjdCwgaWQudmFyPWMoJ0RhdGUnLCAnU3RhdGlvbk5hbWUnKSkKCmBgYAoKCmBgYHtyLCBmaWcud2lkdGg9MTIsIGZpZy5oZWlnaHQ9OH0KCmdncGxvdChkZkNvdmlkRWZmZWN0TWVsdGVkLCBhZXMoeCA9IERhdGUsIHkgPSB2YWx1ZSwgY29sPXZhcmlhYmxlKSkgKwogIGdlb21fcG9pbnQoY29sb3I9J2JsYWNrJyxzaXplPTEpICArCiAgZ2VvbV9saW5lKCkgICsKICBmYWNldF93cmFwKH4gU3RhdGlvbk5hbWUsIG5jb2wgPSAzLCBzY2FsZXM9ImZyZWUiKSArCiAgbGFicyh0aXRsZSA9ICJDb21wYXJhY2nDs24gZGVsIMOtbmRpY2UgZGUgcnVpZG8gTEFlcSBkaWFyaW8gcGFyYSAyMDIwIHkgbGEgbWVkaWEgZGUgbG9zIGHDsW9zIDIwMTUtMjAxOSIsCiAgICAgICBzdWJ0aXRsZSA9ICJTdWJ0w610dWxvIiwgeSA9ICJMQWVxIikgKwogdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJ0b3AiKQpgYGAKCg==